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ABSTRACT 

High resolution numerical simulations of thermal convection in a rap- 
idly rotating channel with gravity perpendicular to the rotation vector are 
described. The convecting columns are subject to a /3-effect resulting from 
cross-channel topographic vortex stretching. The symmetries of the prob- 
lem allow many invariant wavenumber sets, and this property is associated 
with the existence of stable multiple-equilibria at modest supercritical- 
ity. The transition to chaotic behavior involves the production of inter- 
mittent unstable orbits off a two-torus in energy space. At very high Ray- 
leigh number (of order 10‘ to 10 7 ) the motion can be turbulent, depending on 
the size of 0 . However, the turbulence is usually characterized by an 
almost-periodic formation of patches of small scale convection that cause 
regular pulsations in the accompanying strong zonal jets. The processes 
maintaining these flows may be related to those responsible for the zonal 
currents on Jupiter and for cyclic variability on the Sun. 
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0 . OVERVIEW 

The Geophysical Fluid Flow Cell (GFFC) experiment is scheduled to fly 
on board the USML-2 mission in September 1995. In support of this flight, 
and to generate ideas for experiments to be conducted, we have carried out 
an extensive modeling study of thermal convection confined to an equato- 
rial annulus. This is a reasonable model for the nonlinear dynamics of 
quasi-two-dimensional "banana” cells, i.e. the convection modes observed on 
the SL-3 flight of the GFFC. Of great interest are the questions of how the 
banana cell states become chaotic, whether or not mean zonal flows can be 
generated, and what the nature of the turbulent states are. To this end we 
have constructed models that are relatively simple to integrate (as 
opposed to GCM type simulations that are out of the question for the GFFC 
parameter ranges and timescales) . This report describes the major modeling 
results, and focuses on interesting discovery of exotic pulsating turbulent 
states that will be sought during the USML-2 mission. 

1. INTRODUCTION 

One of the outstanding problems in geophysical fluid dynamics is the 
origin of zonal winds on the giant planets and differential rotation on the 
Sun. Busse (1976, 1983) proposed that the Taylor-Proudman constraint asso- 
ciated with a strong basic rotation would cause the interior motions of 
rapidly rotating planetary atmospheres to be invariant along the rotation 
axis. In Boussinesq gas or liquid models the resulting two dimensional con- 
verting columns were found to be subject to a number of secondary instabil- 
ities (Or and Busse, 1987, hereafter OB87, for free boundaries; Schnaubelt 
and Busse, 1991, Schnaubelt, 1992, for rigid boundaries). One is a mean 
flow instability in which converting columns, tilting in the zonal and 
radial directions, produce Reynolds stresses ttiat generate a zonal shear 
flow. This zonal shearing motion is capable of further tilting the cells, 
leading, in principle, to large mean flows as the dissipation is reduced. 

A crucial effect in Busse' s model is the stretching, by the curvature 
of the planetary atmosphere, of radially displaced converting columns 
aligned with the rotation axis. The standard model configuration is an 
equatorial annulus or zonally periodic channel with sloping ends and grav- 
ity perpendicular to the basic rotation vector (figure 1) . The slopes gen- 
erate a topographic /3-effect which provides a mechanism for travelling 
Rossby waves, as well as a tendency for chaotic behavior at modest Rayleigh 
number. In this paper we study the spontaneous generation of strong zonal 
jets at high Rayleigh number. If 0 is either too small or too large, the 
motions are either lack such jets, or are stable altogether. We call the 
intermediate regime, with dominant zonal currents and various forms of 
turbulence, /3 - convection. Although the geometry and vertical structure 
of the Boussinesq model is highly simplified, the results are relevant to 
scaling arguments of Ingersoll and Pollard (1982) for deep convection on 
the giant planets. These order of magnitude estimates rely, in part, on the 
low Rayleigh number computations of Lipps (1971) for ordinary 2-D convec- 
tion in an imposed shear flow. As the presence of vortex stretching associ- 
ated with the /3-effect can have a dramatic influence on the structure of 
both low and high Rayleigh number situations, it is of interest to charac- 
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terize the possible states through a sequence of numerical experiments. 

Thompson (1970) proposed a type of mean flow instability as an explana- 
tion of the 4-day retrograde rotation of the Venusian atmosphere. His 
model was criticized by Lilly (1971) who argued .that on Venus there is no 
constraint to guarantee pure transverse converting rolls that have the 
necessary strong interactions with a zonal shear flow. In the Busse sce- 
nario such a constraint is provided by the basic rotation. Laboratory 
experiments on Boussinesq- liquid spherical shells using outward centrifu- 
gal buoyancy (Carrigan and Busse, 1983) and electrostatic radial buoyancy 
(Hart et. al, 1986) confirmed the validity of this idea over a fairly wide 
range of Rayleigh and Taylor numbers. However, the latter study showed 
that increasing the Rayleigh number at fixed Taylor number invariably leads 
to a breakdown of the columnar "banana cell" convection as the Taylor- 
Proudman constraint is modified by buoyancy forces parallel to the basic 
rotation vector. The breakdown first occurs at high latitudes, then gradu- 
ally erodes down into the equatorial banana cell region as the Rayleigh 
number increases. Thus one question is whether or not there is a parameter 
range for /3— convection with very strong mean zonal flows and self-consis- 
tency in the sense of low local Rossby number and negligible thermal winds. 

There have been a number of studies of the evolution of columnar /3-con- 
vection with increasing Rayleigh number, addressing the question of the 
transition from steadily propagating thermal Rossby waves, which arise at 
the linear stability boundary for the onset of convection, to more compli- 
cated spatio-temporal forms. For application to planetary atmospheres and 
the solar rotation problem, a model should produce a relatively large zonal 
flow in which the columnar convective eddies do not dominate the total 
stream fie Ids . For example, the data from Voyager (Ingersoll, et al., 1981) 
show fluctuating velocities at cloud levels on Jupiter which are consider- 
ably smaller than the zonal-mean jet speed at most latitudes. However, all 
the weakly nonlinear models of /3-convection constructed to date have eddy 
kinetic energies either far exceed the zonal means. Highly truncated 
models have been pushed to large Rayleigh number where rough equivalence 
can be found, but the use of low truncation at high super criticality is sus- 
pect. In both cases strong regular-wave patterns at modest wavenumber are 
found that are not easily associated with features imbedded in the clouds 
on the giant planets, nor in photometry or helioseismology of the solar 
atmosphere. All these objects nonetheless boast a strong latitudinally 
varying differential zonal rotation (Howard and Harvey, 1970, Ingersoll et 
al., 1981). 

Lin, Busse and Ghil (1989), hereafter LGB89, used a perturbation method 
to study mode mixing and transition, ostensibly near the neutral curve for 
linear instability. Their model predicted two "multiple" states with the 
same waveset symmetry but differing amounts of zonal cell tilt with radius. 
Transition to chaotic time dependence was via a pe riod dou bling cascade at 
near— critical Rayleigh number. Lin (1990) used a truncated Fourier expan- 
sion to suggest that the convecting column system might be subject to a 
sideband "double column" instability in which rows of like-sign vortices 
would align permanently on the same side of the model channel, leading to 
large mean flows at modest Rayleigh numbers. A period doubling cascade was 
predicted at low super criticality in this latter work as well. It is not 
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clear that these results are robust with respect to model truncation. The 
well-known failure of severely truncated models of two-dimensional con- 
vection and baroclinic instability to predict the outcome of highly 
resolved simulations at even modest supercriticality has motivated us to 
revisit the thermal Rossby wave problem. In our high resolution calcula- 
tions we find significant differences in the transition scenario compared 
to the low-order models, as well as some unusual turbulent states. There 
are interna lly-consistent turbulent regimes with large zonal jets that may 
either be quasi-steady or almost-periodic, depending on 0 . 

The rest of this report is organized as follows. The basic model equa- 
tions and numerical methods are summarized in section 2. Section 3 
describes our results on multiple equilibria and the transition to chaos. 
Section 4 presents data from a number of runs at large Rayleigh number with 
various values of the 0 parameter. Our main conclusions are summarized in 
section 5. 

2 . MODEL EQUATIONS AND METHOD 

Our central motivation is to shed light on the nonlinear dynamics of 
/3-convection in the simplest geometry that has been universally used in the 
previous modeling studies mentioned above. The basic model is derived for 
flow of a Boussinesq liquid in a rectilinear channel whose cross-section is 
shown in figure 1. The governing equations are obtained by perturbing off a 
state of qeostrophic balance (Busse and Or, 1986). The Rossby number for 
the flow, U/20D or, as an even stronger condition, w*/20, is assumed small. 
Here U is the rms flow speed in the x-y plane, 0 is the basic rotation rate, 
D the channel width, and «* is the local vertical vorticity. The motions 
are then essentially two-dimensional in the x-y plane and controlled by the 
quasi-geostrophic vorticity equation (c.f. Pedlosky, 1987), with the one 
addition that zonal (x) variations of the relatively weak y-directed buoy- 
ancy forces can generate "vertical" vorticity (i.e. vorticity aligned with 
the rotation axis) . In order to maintain horizontal non-divergence of the 
lowest order velocity field in the x-y plane, the component of buoyancy in 
the y-direction is assumed small compared with the Coriolis deflection of 
the zonal u velocity. The thermal winds associated with horizontal temper- 
ature gradients and the component of gravity along the axis of rotation 
that would be present in a spherical atmosphere are neglected. In order to 
be self-consistent with lowest order geostrophy, the nonlinear vortex 
stretching term in the vertical vorticity equation is neglected with res- 
pect to planetary vortex stretching. These scaling assumptions are re- 
examined a' posteriori in section 5. Finally, while topographic slopes are 
used to generate a /3-effect, all the boundaries are assumed to be free- 
slip. If the sloping surfaces are almost flat, the lowest order fields do 
not require Ekman layers to adjust the normal shear of the horizontal velo- 
cities at these surfaces to zero. Thus there will be no Ekman suction damp- 
ing, and lateral diffusion of heat and momentum is the only dissipative 
mechanism included. 

When lengths are scaled by D, velocities by rD” 1 , temperature relative 
to a background conductive state by i-r'AT, and time by D 2 v~‘ , the governing 
vorticity equation and the heat equation become (following LBG89): 
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JH + J »'“> - » M = - Ra “ + ^ 

P r {®«(W)} = fT-“, (2) 


1 — 

with 

J 

w — V 2 \p , 

(3) 


u-- 

ay 

(4) 

and 

v= ft . 

dx 

(5) 

The parameters are defined by: 



Prandtl number P r * L , 

K 


(6) 

Rayleigh number R a « f 


(7) 

0 - parameter 0 m AtanXaiflS . 

JuV 


(8) 


Here a is the coefficient of thermal expansion, v the kinematic viscosity, 
and k the thermal diffusivity. Other parameters are defined in figure 1. 
This paper studies flow regimes in the R a - 0 plane at a fixed value of P r = 
1 . 0 . 


The 0 parameter as given in (8) is related to the topographic slopes in 
figure 1, which in turn follow from extracting a slot from a full sphere 
and assuming that the motions extend along the axis of rotation all the way 
through the interior. This was the canonical view proposed by Busse (1983). 
However it is important to observe that eqn. (8) is identical in magnitude 
to the 0 parameter obtained for shallow to moderately deep convection 
between spherical shells separated by a distance H, measured along the axis 
of rotation, when H is small compared to a full sphere radius R s . Under the 
local slope approximation for a slot of width D intersecting the surface of 
a sphere of radius R s at latitude 6, the ratio of the change in depth across 
the slot to the average depth is 

G m £& _ 20D IQ m ItegsM (8 , j 

R s sin J (0)|l - — — ? . . 1 
s '{ RgSin^) J 
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where the last term in brackets is absent for columnar convection extending 
all the way through. Since the column length L (or H) only appears in 0, the 
results for a given 0 apply equally to shallow or moderately deep as well 
as to full interior-crossing Boussinesq convection. 

The above system is solved numerically by a pseudospectra 1 algorithm. 
The expansion functions are chosen such that stress-free, impermeable, and 
infinitely conducting boundary conditions are satisfied on the vertical 


walls at y = 0 and y = 1. We 

write: 



00 

00 



I 

Z 

sin (my) exp(imx) , 

(9) 

m=0 

n=l 



00 

00 



*- I 

V T 
A, A n#n 

sin (my) exp(/mx) , 

(10) 

m=0 

n=l 




Then, the boundary conditions mentioned above, 

^=2M=T=0, at y = 0,1, 

and periodic zonal boundary conditions are automatically satisfied. The 
above representation assumes that the zonal periodicity length is 2x times 
the channel width. This particular aspect ratio is maintained throughout 
the paper. The solutions are affected if the aspect is significantly 
smaller so that the lowest wavenumber corresponds to tall thin cells. 
Linear stability analyses of finite-amplitude /3-convection in an infinite 
channel (Schnaubelt, 1992) predict sideband instabilities with maximum 
growth rates for wavenumbers corresponding to m = 1 in a 2ir long channel. 
Thus our choice of this particular aspect ratio is expected to capture the 
essential physics of flows in longer channels while maximizing computa- 
tional efficiency. 

The numerical method is a Fourier spectral collocation scheme based 
upon a truncation of (9) and (10). The variables \p and T of configuration 
(physical) space are mapped back and forth into the coefficients ^ m/n and 
T m n of spectral space using fast Fourier transforms. The time integration 
of ’equations (1) and (2) is performed semi- implicitly. An implicit Crank- 
Nicolson formulation (Smith, 1965) is used for the linear terms, whereas an 
explicit three-level Adams-Bashforth scheme (Richtmyer and Morton, 1967) 
is adopted for the nonlinear terms. The Crank-Nicolson step is uncondi- 
tionally stable and therefore the timestep for the linear parts is 
restricted solely by requirements of accuracy. However, the explicit 
scheme for the nonlinear parts requires a timestep that satisfies the Cour- 
ant-Friedrichs-Lewy condition for numerical stability. 

Products of the variables are evaluated in configuration space to 
avoid computationally expensive convolution sums of the phase space vari- 
ables. However, all other calculations are carried out in phase space, 
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where derivatives reduce to simple multiplications by powers of the waven- 
umbers. Aliasing errors are removed using the so-called 2/3 dealiasing 
rule. The pseudospectra 1 code allows initial conditions to be random fluc- 
tuations in any required wavenumber set, or the results of a previous simu- 
lation. It was thoroughly tested against linear analysis, simple nonlinear 
examples, and (by setting 0 = 0) against previous two-dimensional Rayleigh- 
Benard results (e.g. Moore and Weiss, 1973). The results presented below 
are believed to be convergent on the basis of obtaining quantitatively sim- 
ilar solutions at double the spatial resolution. The higher Rayleigh 
number runs used 256 (x) and 65 (y) modes. A typical run of 200,000 times- 
teps takes about a day on an IBM RS6000-320 workstation, or a couple of 
hours on a CRAY-YMP. 

3 . MULTIPLE STATES AND THE TRANSITION TO CHAOS 

We start by summarizing a large number of computational nans at moder- 
ately supercritical Rayleigh numbers designed to illustrate the types of 
secondary instability involved in the transition to chaos. We also show 
that in this intermediate R a region of parameter space between linear in- 
stability and full turbulence there are at least two locations where line- 
arly stable multiple solutions are found. All the results reported in this 
section are for P r = 1.0 and 0 = 2800, to permit comparisons with previously 
published works on this problem that have concentrated on these particular 
values. 

As a point of reference, figure 2 illustrates the behavior of the pri- 
mary instability. The linear versions of (1) - (2) have solutions with $ and 
T proportional to exp(st + /ax) sin (my) . The growth rates for P r = 1 are 
then given by the quadratic 


k?} " + k4 + i0a = 0 ' (11) 

with k 2 ■ n 2 ^ + a 7 . From this one can easily determine the critical curve 
where Re(s) crosses the imaginary axis. This occurs at a frequency 


s 2 + s( 2k 2 + 


« 


= zM 
2k 2 


and a critical Rayleigh number 


R ac 



JL 

4k 2 


( 12 ) 


(13) 


Figure 2a shows that the critical Rayleigh number and wavenumber both in- 
crease with 0 1 (as 0*h and 0'h respectively, for 0 large). At supercritical 
conditions the wavenumber of maximum growth is smaller than the critical 
wavenumber (e.g. compare figures 2a and 2b for 0 = 20000). In our calcula- 
tions the x-wavenumber a = m is an integer quantized by the zonal periodi- 
city. For P r = 1.0 and 0 = 2800 the critical values are R ac = 30,830 at a c = 
m = 9 and n = 1. The critical Rayleigh number for n = 2, m = 9 is slightly 
bigger (36,254), but as n increases beyond 2 the critical Rayleigh number 
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increases rapidly (for n = 3 it is about 67,000) and the associated critical 
wavenumber first decreases then increases. Over the whole parameter plane 
0 < R a < 10 s , 0 < 0 < 10« the n = 1 mode has the highest linear growth rate, 
although as 0 becomes very large a wider and wider range of n' s have roughly 
the same growth rate. This leads to the important question of finite-amp- 
litude meridional (i.e. cross-channel) scale selection. 

Nonlinear solutions of (l)-(2) can be found with certain invariant wav- 
enumber subsets based on zonal and cross-channel periodicity, and on a 
form-preserving shift-ref lect symmetry. It is obvious that any solution 
with initial zonal wavenumbers confined to m = 0 (the zonal component of 
the flow), m = M (the "fundamental wave"), and all zonal harmonics of M, 
will not scatter energy into any other zonal wavenumbers. Solutions based 
on dividing the zonal periodicity length by some integer M thus preserve 
the initial sparse (for M large) spectral occupation matrix, though they 
may be unstable to perturbations in longer wavelengths. For example, OB87 
study linear perturbations to finite amplitude solutions in an infinite 
channel with various values of M around 8, and find supercritical sideband 
instabilities with long wave excitation in wavenumbers smaller than one. 
Similar comments can be made with respect to the cross-channel periodicity 
(for these boundary conditions). However, we do not find stable finite- 
amplitude solutions with a fundamental cross-channel periodicity N > 1, 
though they certainly exist. 

Invariant wavenumber sets with only even values of n + m in (9) and (10) 
are also possible (OB87). These are related to motions in the channel that 
are preserved under a reflection about the mid-line y = 0.5, and a transla- 
tion in x by one half the fundamental M“‘ periodicity length. In the nonli- 
near solutions studied here we label various end-states based on their fun- 
damental zonal periodicity M, and whether they have the shift-ref lect sym- 
metry or not. Asymmetric A-states have both even and odd values of n + m, 
while S-states have n + m even only. In addition, it is useful to note that 
some states, while filling low wavenumbers, still have peak energies at a 
fairly large value of m. Figure 3 shows spectral occupation diagrams for 
several different states and illustrates our nomenclature. Finally, every 
solution has a companion reflected about the midplane under y -► 1 - y 
alone, and/ or shifted by an arbitrary phase in x. The simplified geometry 
with uniformly sloping ends does not provide enough asymmetry to determine 
the sign of the zonal flow at y * 0, for example, and this, as well as the 
eddy positions in x, are determined by initial conditions. The above symme- 
tries are often found in such 2-D channel models. Some examples are ordi- 
nary sheared convection, Howard and Krishnamurti (1986), baroclinic insta- 
bility on the f-plane, Cattaneo and Hart (1990), and thermosolutal convec- 
tion, Moore et al. (1991), among others. Because of the 0-effect our prob- 
lem does not share the left-right x -* -x reflection symmetry that these 
latter systems have. In addition it is useful to note that the shift-ref - 
lect symmetry S-states will be lost in more general geometries such as cyl- 
inders, spherical shells, bodies with non-conical ends, etc. As shown 
below, S-modes participate in the low Rayleigh number bifurcations, but are 
unstable at suitably high forcing. 

Figure 4 summarizes the results of about 100 runs with P r = 1.0 and 0 = 



- 10 - 


2800. At slightly supercritical values of R ? , the system equilibrates to a 
type 9S symmetric state, as expected from linear theory. At a value of R a 
slightly smaller than 35,000 this 9S state becomes unstable to the 9A mode, 
consistent with the finite amplitude stability calculations. Or and Busse 
(1987) call this symmetry breaking instability a "mean flow instability", 
because the interaction of zonally phase-shifted even and odd cross-chan- 
nel inodes with the same zonal wavenumber can generate an m - 0 zonal flow 
component. LBG89 give examples of multiple equilibria following the bifur- 
cation to an asymmetric state. Only one asymmetric type 9A solution was 
found in the present computational exercise. The phase shift between the 
first and second cross-stream mode of the fundamental zonal wave is 130 
degrees at R a = 35,000. This corresponds to "Solution I" of LBG89, that has 
a phase shift of "x = 120°". We attempted to start a 9A solution at Rayleigh 
numbers near 35,500 with an initial phase shift of 15° and relative ampli- 
tudes given by LGB89 for their "solution n". However the system typically 
evolved to a singly periodic mean flow state with x — 130°. This corres- 
ponds to the vacillating solution I of LGB89. Based on our limited and 
potentially protocol-sensitive study of this issue, we suspect that there 
is only one asymmetric solution at near critical conditions. 

The asymmetric 9A state bifurcates to an amplitude-periodic 8A state 
at R a = 37,500. OB87 showed that the mean flow secondary state (i.e. 9A) 
would almost immediately (upon increasing R a ) suffer a sideband instabil- 
ity, but they did not find any tertiary steady convection states. Thus we 
believe the 8A limit cycle to be a simple extension of the 9A per iodic- state 
similar to the so-called vacillatory convection solution calculated by 
OB87 for R a near 40,000 at low resolution (triangular T5 truncation; 
retaining only modes with n + m < 6). Because they did not consider finite 
aspect ratio channels, the comparison cannot be made more exact. The sim- 
ilarities suggest that periodic amplitude modulation of the convection, as 
seen in the limited resolution model of LGB89 and the O(10 2 ) degree of free- 
dom calculations of Schnaubelt (1992), is a robust feature of this problem. 

However, there is simultaneously a regime of stable steady convection 
following the instability of the mean-flow 9A mode, that has, however, the 
1A(7) symmetry. This nonlinear fixed-amplitude tertiary travelling wave 
state is a result of the sideband instability. The 8A and 1A branches coex- 
ist for a substantial range of R a . Initial value problems with small 1A 
perturbations to the 8A state indicate that it is stable up to R a a 60,000 
which is close to the onset-value for totally desymmetrized chaos. Each 
integration at increasing R a along the 1A branch is essentially subject to 
a 1A perturbation (i.e. the previous solution), and so is at least linearly 
stable until a transition is made to another type of time dependence or 
waveset symmetry. 

Previous studies have cited evidence for period-doubling cascades to 
chaos. At high resolution, completed cascades of this type are not found. 
The 8A state undergoes a single doubling at R a = 57,500, but then a further 
bifurcation to a torus occurs in a narrow window centered on 75,500, before 
chaotic motion ensues for R a > 76,000. Figure 5 illustrates this transi- 
tion. Note that the single period-doubling occurs at much higher values of 
R a than that found from highly truncated models with period doubling (R a = 
36,000, LGB89 , Lin 90), indicating that previous studies used insufficient 
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resolution to explore the actual transition to chaos. OB87 also found 
period doubling in this range, but used somewhat different boundary condi- 
tions having free-stress eddies, but periodic zonal flows. It is not known 
whether the differences between our calculations and theirs are due to 
boundary conditions or to truncation. In our 8A transition, the initial 
high frequency mean-flow vacillation is eventually replaced by a low-fre- 
quency motion that has its origin in the toroidal bifurcation to modulated 
vacillating convection. Figure 6a shows that the limit cycle behavior is 
related to a periodicity in the tilt of the convection cells that drives an 
oscillating zonal flow. This is largely a result of a periodic phase shift 
in x of the even and odd cross-stream fundamental modes, giving rise to a 
sort of interference vacillation. At larger R a there is a near ly-per iodic 
pairing of like sign vortices on opposite sides of the channel (e.g. figure 
6b panels 5 and 13), due to the relative drift of the row of vortices at the 
top of the channel with respect to those at the bottom. This is reminiscent 
of the truncated solutions of Lin (1990) where rows of n = 2 vortices pro- 
pagate relative to each other along opposite walls at modest R a . Lin pred- 
icted that the like-sign vortices would permanently position themselves 
near the opposing boundaries as R a is further increased, in a process that 
was called a "double column" instability. The solution in figure 6b might 
more appropriately be called an n = 2 differentially-dr if ting-column (DDC) 
state. Similar vacillations were predicted by Schnaubelt (1992) in situa- 
tions with curved ends, giving a y-varying 0 parameter. He found oscilla- 
tory instabilities of steady convection with the DDC signature. The DDC 
vacillation appearing here at high R a does not require sidebands or curved 
boundaries for its existence. The intermittent excursions off the torus 
(figure 5c) are suggestive of a tangent bifurcation. Physically they seem 
to be related to almost periodic pulsations of the zonal flow - convection 
system, which, because of small scale high-frequency fluctuations, can 
occur irregularly in time. 

For completeness we note that the IS state undergoes a similar transi- 
tion to chaos, though at lower R a . Referring to figure 4, the steady 9S 
mode first suffers a sideband instability in which the new the spectrum is 
dominated by wavenumber m = 9. With increasing R^ the spectrum flattens 
out and the flow makes a transition to aperiodicity at a low Rayleigh 
number, R a = 43,750. Figure 7 emphasizes that there is no period doubling, 
and the frequency separation in the quasi-periodic states is not large. 
Because they are unstable, both the IS and the 8A transition branches are 
not physically significant. 

For R a > ==62,500 the only stable state is 1A(6). The chaotic branch 
overlaps with a steady solution in a narrow window between 59,500 and 
62,500. The steady solution is dominated by wavenumber 7 and has much 
weaker sidebands than the wavenumber 6 dominated chaotic solution at the 
same R a . Figure 8 shows typical streamfunctions and zonal wavenumber 
spectra for these modes. The spatial spectra of the 1A chaotic branch, 
though dominated by m = 6 near R a = 60,000, rapidly become filled in at low 
wavenumbers as R a increases. We could not track the 1A chaotic state back 
in Rayleigh number to find its origin. At R a = 59,500 it loses stability to 
a 2A periodic state, that in turn becomes unstable to 8A at slightly lower 

R a • 
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4. HIGH RAYLEIGH NUMBER FLOWS 

The zonal (m = 0) component of the flow over the regime diagram of 
figure 4 does not exceed the eddy velocities and is generally rather small. 
It is of interest to explore larger Rayleigh numbers to identify parameter 
values where one might get substantial zonal jets. We begin with a cut in 0 
at R a = 10* , P r = 1.0 . Some typical streamf unction vs. time images are 
shown in figure 9. These solutions are obtained from desymmetrized initial 
conditions (type 1A) , and all but panels a) and b) remain asymmetric. 

At 0 - 0 one recovers steady non-chaotic roll convection dominated by 
wavenumber 1. A small (0(10 3 )) value of 0 leads to steady travelling con- 
vection, with waves travelling from left to right in response to the 
deepening of the channel in y. 0 ’ s of this magnitude produce weak mean 
flows with space-averaged mean kinetic energy K much less than the eddy 
kinetic energy K' . Near 0 = 6000 there is a narrow window of hysteresis 
between travelling convection with almost-steady amplitude and weak mean 
flows connected to the 0 = 0 ordinary convection branch, and a chaotic 
strong- jet jet solution that is connected to the turbulent /3-convection 
branch. The hysteresis spans a 0 range of several hundred (see also fig. 
14). As 0 is further increased the dominant wavenumber of the convection 
grows, but it remains substantially less than that expected from linear 
theory. Figure 10 shows the relation between nonlinear wave selection and 
the fastest linearly growing wave for the same conditions. The difference 
is biggest at low 0 where the flows are most supercritical. The largest 
zonal flows are found on the low 0 side of the chaotic domain, and in this 
region the zonal circulation is quasi-steady. When 0 is increased almost- 
periodic vacillations of the zonal shear flow and the heat flux across the 
channel are observed (figure 11). 

Figure 12 shows total stream and temperature fields (less the 
conduction profile), which are typical of flows at extreme R a . The associ- 
ated time series of zonal energy, pointwise zonal velocity and heat flux 
are shown in figure 13 The zonal flow is built up rapidly during a p erio d 
of violent and small scale (because 0 is large) convection. When the zonal 
shear reaches its maximum the convection is suppressed and the zonal jet 
relaxes slowly back to low values and the process starts over. The rela- 
tively slight aperiodicity of this almost-periodic large scale system 
appears due to the irregularity and patchiness of the small scale turbulent 
convection that arises when the shear is small. The cycling is consistent 
with the idea of mean flow instability and rapid zonal shear growth as the 
convection turns on, followed by the relatively quick suppression of the 
convective plumes when the shear becomes too large (we recall that linear 
stability results show stabilization of transverse roll convection by a 
strong vertical shear) . An important point is that the zonal shear relaxa- 
tion time is roughly independent of 0 and R a (see figures 11 and 13 and 
table 1), and therefore this turbulent flow vacillation occurs on the vis- 
cous D 2 /*' timescale. Once the convection shuts down, the zonal jets decay 
on this long time scale. The mechanism for these cycles appears to be dif- 
ferent from the mode-mixing interference vacillations at low R a (i.e. fig. 
6b). 
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TABLE 1 

CHARACTERISTICS OF 0-CONVECTION 


0 

R a 

Nu-1 

K 

K' 

AU 

K'K 

r 

0.0 

5x10* 

15.1 

0 

2x10* 

— 

_ 

_ 

4000 

It 

7.1 

2.8x10* 

8.7x10* 

22 

6.3x10* 

— 

4500 

VI 

4.5 

1.8x10* 

5.6x10* 

76 

2.3x10* 

— 

4500 

VI 

4.3 

1.3x10* 

9.9x10* 

702 

8.4x10* 

.12 

5000 

IV 

1.9 

1.1x10* 

5.4x10* 

654 

4.2x10* 

.12 

1.0X10* 

IV 

1.4 

1.5x10* 

4.0x10* 

163 

1.4x10* 

.12 

1.5X10 4 

VI 

1.0 

2.3x10* 

6.3x10* 

64 

2.8x10* 

— 

0.0 

lxlO 6 

19.1 

0 

4.9x10* 

— 

- 

- 

4000 

VI 

16.4 

3.5x10* 

4.2x10* 

17 

7.6x10* 

— 

6000 

VI 

9.1 

5.6x10* 

2.1x10* 

25 

2.1x10* 

— 

6000 

VI 

5.4 

6.2x10* 

1.6x10* 

1526 

3.2x10* 

- 

8000 

VI 

1.8 

1.9x10* 

6.9x10* 

880 

1.1x10* 

.12 

1.5x10* 

IV 

.77 

5.1x10* 

1.5x10* 

479 

2.7x10* 

.08 

2.0x10* 

IV 

.51 

2.3x10* 

7.1x10’ 

326 

1.2x10* 

.08 

2.8X10* 

IV 

- 

2.0x10* 

5.4x10* 

12 

4.0x10* 

.06 

2 .8X10* 

5x10* 

3.67 

9.6x10* 

3.4x10* 

224 

2.7x10* 

.08 

3.5X10* 

vv 

1.59 

5.6x10* 

4.1x10* 

1588 

2.6x10* 

.07 

5.0X10* 

VI 

1.05 

2.4x10* 

4.4x10* 

1056 

1.5x10* 

.12 

7.5x10* 

VI 

.47 

7.2x10* 

1.6x10* 

577 

4.1x10* 

.09 

1x10* 

VI 

.26 

1.4x10* 

6.3x10’ 

252 

7.3x10* 

— 

4.0x10* 

lxlO 7 

3.7 

1.7x10* 

7.6x10* 

315 

9.4x10* 

.07 

4.5x10* 

VI 

2.3 

1.8x10* 

1.7x10* 

2839 

1.1x10* 

.10 

5.0x10* 

IV 

1.8 

1.5x10* 

1.8x10* 

2581 

6.4x10* 

.10 

7.5x10* 

IV 

.93 

6.9x10* 

1.7x10* 

1803 

3.1x10* 

.08 

1.0x10* 

II 

.73 

3.3x10* 

5.4x10* 

1264 

1.6x10* 

.07 

1.5X10* 

IV 

.33 

6.3x10* 

1.5x10* 

552 

3.1x10* 

— 


Nu - 1 mx and time average of w'T' at y =0.5. 

K ■ global and time average of zonal kinetic energy (m -0.) 

K' m global and time average of eddy kinetic energy (m *0). 

A U mtime average of zonal cross- channel velocity difference. 

K'K * space- time average of eddy to zonal kinetic energy production rate. 
t ■ dominant period of zonal kinetic energy fluctuation. 

(shorter run times lead to uncertainties in some lines of the table) 


Figures 14, 15 and table 1 summarize the situation. As R a increases a 
window in 0 opens up where large zonal jets can be found. Quasi-steady 
zonal currents occur on the low-0 side of the window, while almost periodic 
jets occur in the middle. At fixed R a , increasing 0 causes a drop in the 
convective flux N u - 1, along with a concurrent decrease in the eddy kinetic 
energy. The zonal kinetic energy first rises then falls, and the largest 
values of the shear AU/D seem to occur when the eddy and zonal kinetic 
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energies are roughly the same, especially at high R a . However, there are 
many states with K substantially greater than K' and these occur as 0 moves 
towards, but is not too close to, the linear stability boundary. Further 
discussion follows in section 5. We have not taken the extremely small 
steps in 0 needed to explore the high Rayleigh number 0 -» 0]^ transition in 
detail, as simulations of these low-friction motions require high spatial 
resolution and small timesteps. 


5. DISCUSSION AND CONCLUSIONS 

We have studied two-dimensional thermal convection aligned with the 
axis of basic rotation in a 0-plane channel of aspect ratio 2 t having free 
slip sidewalls. The computations had unit Prandtl number. The main conc- 
lusions are: 

i) The fundamental single-wave state, which arises from the linear in- 
stability of the conductive thermal profile for 0 less than 


^1 in " 





f 


(14) 


equilibrates at wavenumber 9 (18 cells in the zonal direction for 0 = 2800). 
As the Rayleigh number is increased at fixed 0, this finite-amplitude 
steady convection suffers successive secondary instabilities to an asymme- 
tric (mean flow) steady state, and then to a vacillatory mean flow state at 
a slightly lower wavenumber (8). This follows the scenario of Or and Busse 
(1987), who first performed linear stability analyses of the finite-ampli- 
tude supercritical steady states. 

ii) The vacillatory motion, including only wavenumber 8 and its harmon- 
ics, becomes chaotic at about twice the critical Rayleigh number via a 
quasi-periodic 2-torus scenario. The two-frequencies seem to result from a 
tilted cell interaction and a phase-winding instability of the first two 
cross-stream modes. Earlier studies (Lin et al. 1989, Lin 1990) that pred- 
ict a transition to chaos by period-doubling at about 20 percent supercri- 
ticality may be inaccurate due to truncation. Another possible cause for 
the difference may be related to the finite aspect ratio geometry with its 
associated wavenumber discretization. 

iii) Only 1 steady asymmetric form of convection was found at slightly 
supercritical Rayleigh number. The zonal phase relations between the first 
two cross-stream wavy modes suggest that this state is consistent with the 
"Solution I" weakly-nonlinear perturbation results of Lin et al. 1989. The 
absence of their "Solution H" may be a result of our numerical protocol. 

iv) For 40,000 < R a < 60,000, and 0 — 2800, several linearly stable 
finite amplitude states coexist. These include a periodic wavenumber-8 
(plus harmonics) state, a periodic wavenumber-2 (plus harmonics) state, and 
a steady asymmetric state with all wavenumbers excited. This last mode is 
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the roost important because of its association with the high Rayleigh number 
fully desymmetrized turbulent flows. 

v) For Rayleigh numbers up to 10 7 the most vigorous zonal jets induced 
by the convection can be found in parameter space along the long-dash line 
of figure 15. This is approximately given by 


0v - 



(15) 


For 0 < 0 V the amplitude of the zonal velocity falls rapidly, and is zero 
when 0 = 0. In this region the motions are weakly chaotic in a band with 0 
slightly less than 0 W , becoming steady in amplitude for 0 less than about 
0.5 0 V . This latter value is a crude estimate based on a coarse sampling in 
this area of parameter space. Along the 0 V curve the dominant zonal waven- 
umber of the convection is one or two, and the zonal flows are quasi-steady 
with small but irregular fluctuations associated with the high-frequency 
structure in the convection. 

For 0 > 0 V the dominant wavenumber of the convection increases and the 
zonal flows begin to pulsate on a diffusive timescale. The small scale 
convection is highly turbulent, especially at our extreme value of R a = 10 7 . 
However, the oscillation of the zonally averaged quantities is nearly regu- 
lar. 


Along a curve very crudely given by 

0m * *013 R a (16) 

the ratio of the time-averaged zonal kinetic energy to the time-averaged 
eddy kinetic energy is a maximum, with values of about 2.5 at R a = 5xl0 5 and 
over 6 at R a = 10 7 . It is doubtful that intersects the linear stability 
boundary (14) though it is probable that the turbulent pulsating states 
occur closer to 0Hn as R a is increased beyond 10 7 . As 0 is increased above 
0m at fixed R a the convection shifts to still higher wavenumber and eventu- 
ally dies out as the neutral curve (14) is crossed. 

If R a is increased at fixed 0 (as in section 3), curve (15) is eventually 
crossed and the chaotic convection at smaller R a becomes regular again. 
This process is distinct from the relaminarization found in low order 
models that give steady convection following an inverse period-doubling 
cascade at fixed wavenumber as R a crosses about 40,000 at 0 = 2800 (LGB89) . 
Here the chaotic high-mean-flow states at 0 = 2800 die out at R a » 300,000 
after the flow has made many transitions towards wavenumber 1. This loss 
of turbulence at high R a is an artifact of the 2-D model, which will break 
down if R a is made too large while 0 is held fixed. Such behavior was seen 
in the spherical laboratory experiments of Hart et al. 1986 that show 
strong 3-dimensional convection in this limit. 
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vi) The model appears to be reasonably self-consistent over most of 
its parameter range. The basic expansion requires that the Rossby number 
be small. The root mean square value of Rq over space and time is Rq = 
Go>rins /0 ' where w is the local nondimensional vorticity in the solutions. 
As shown in §2 the parameter G is geometrically related to the actual topo- 
graphic slope in figure 1 by 

G = 2 t a a falP = 2 SB &Jm * o(l) (17) 
L L 

for a planetary atmosphere at subtropical latitudes_ 6 (see also eqn. 8 '). 
Setting G = 1, at R a = 10 * we find, for example, that Rq falls below one for 
0 > 6000 and below one tenth for 0 > 10 4 . 


Horizontal non-divergence of the lowest order velocity field requires 

that 


1 » 


gall: = R s I 

20U* a 0 u * 


(18) 


Using time and space averaged values for temperature T and velocity U, this 
condition is well satisfied by our solutions. For example, at R a = 5xl0‘ 
the RHS of (18) is 0.071 at 0 = 4000, decreasing to 0.021 at 0 = 20,000. 


On a spherical planet thermal winds giving a variation of velocity 
along the axis of rotation will be induced by the component of gravity gp 
parallel to the rotation axis. Combining this component with the slopes 17 
appropriate to a sphere yields an estimate of the fractional change of vel- 
ocity over the depth D be small. This is 

2R a . . 

€= 1#I vt I' (w>. 

since LGgp/gD = 2Dtan(0)/L • cot ( 0 )L/D = 2. Computations of 

(VT.VT) 1 / 2 / (v.v ) 1 / 2 indicate that the e is rarely above 1/2 except at special 
points where the speed is zero. Thus the computed flows appear approxi- 
mately consistent with the model assumptions provided one is in the 0-con- 
vection regime. We expect some quantitative corrections when ageostrophic 
and thermal wind effects are included. At values of 0 <= 0 V the solutions 
may be quite unrealistic for situations with G * 0(1). However, these 
latter motions could be observed in a laboratory centrifugal annulus 
experiment with almost horizontal ends (ij -* 0) so that G becomes very 
small. 


vii) The basic model is too simplistic for a direct quantitative com- 
parison with planetary and stellar atmospheres as it neglects full spheri- 
city, compressibility, as well as other potentially important thermody- 
namic effects. Nonetheless, it is tempting to speculate that Jovian 
dynamics are associated with a turbulent flow having a 0 near to / 3 V , and 
that the convection zone on the Sun is associated with flow near to 0 m . 
This puts Jupiter in a low-wavenumber eddy state associated with a strong 
quasi-steady zonal flow. The Sun is represented by a higher wavenumber 
eddy state with a pulsating zonal flow (e.g. fig. 9c or 9d for Jupiter, fig. 
9e or 12 for the Sun). Magalhaes et al. (1990) have deduced the presence of 
slow moving zonal wavenumber 9-11 eddies in thermal images of Jupiter. 
They suggest that these motions may be rooted many scale heights below the 
cloud deck. On the other hand, surface observations and helioseismology of 
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the interior of the Sun have not as yet exposed any low wavenumber giant 
cells , so we believe that the finer structure high-/? convection may be more 
appropriate for this body. 

After fixing the geometry and the Prandtl number, the model has two 
virtually unknown dimensional parameters, the small-scale eddy viscosity v 
and the super-adiabatic buoyancy forcing gaAT/D exciting the convection. 
Assuming a value of R a = 10 7 and N u = 2, we can obtain an small-eddy diffu- 
sion velocity for Jupiter by requiring that the model yield the observed 
internal convection flux F of about 5 watts/m 2 . The diffusion velocity 


'/ 3 

L _ J qgDF 1 

D { pC p R a N u j 


( 20 ) 


obtained in this way is relatively insensitive to the various parameters. 
Taking Cp = 1.4xl0 4 joules/Kg°K , a = 1/200°K , g = 20 m/sec 2 , p = 2 Kg/m 5 , and 
D = 2xl(r m, we get t>f D = .045 m/sec. Setting G = 0.3 and using the rotation 
rate for Jupiter we get /? = 48,000. The results from table 1 then give a 
dimensional cross-channel zonal velocity difference of about 120 m/sec, 
which is roughly the same size as that observed. The model zonal shear 
fluctuates by about 15 percent on a timescale of 1.3 years. The zonal kin- 
etic energy generation rate, 6.4xl0 6 , is about one third of the eddy kinetic 
energy generation by buoyancy work R a N u = 1.8xl0 7 . The former number tran- 
slates to a dimensional value of 3xl0 -5 m 2 /sec 3 , which is somewhat smaller 
than the values reported from Voyager I of (1 - 4)xl0” 4 m 2 /sec 5 (Ingersoll et 
al. 1981). It may be noted that Sorovsky et al. (1982) have argued that the 
actual correlation, after accounting for sampling bias, may be considerably 
smaller than this latter value. Nonetheless, the model efficiency of 
/3-convection for generating zonal jets is quite high, with a substantial 
fraction of the buoyancy work going in maintenance of the zonal jets. 

Carrying out the same procedure for conditions typical of the outer 
convection zone of the Sun with D equal to a solar radius gives a very large 
viscosity v - 1.3x10’ m 2 /sec. The associated /3 =7500 is below the /3-convec- 
tion regime at R a = 10 7 . If we take a value v = 10* m 2 /sec, which is about a 
third that estimated from surface magnetic feature diffusion and typical 
of values used in large eddy models of compressible convection on the Sian 
(Gilman and Miller, 1986), then 0 = 80,000 and we are in the pulsating tur- 
bulence regime. The diffusion velocity is .14 m/sec and the dimensional 
zonal cross-channel velocity difference is about 300 m/sec, which is the 
observed scale of differential rotation between the equator and high lati- 
tudes. More interesting is the pulsation period, which turns out to be 13 
years for these parameters. Might almost-periodic shear-convection inter- 
actions in a turbulent flow have something to do with the ll year solar 
activity cycle? The demonstration of a dynamical mechanism for producing 
oscillations with a similar timescale is intriguing. However, this latter 
more realistic value of the viscosity, although yielding reasonable 
numbers for the differential rotation rate and the pulsation period, leads 
to a convective flux that is too small by a factor of a thousand. Leaving 
aside the question of flux partitioning between different transport 
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mechanisms and between various represented and non-represented eddy sizes, 
it is thought that a Prandtl number considerably smaller than one is more 
appropriate for the Sun. Low Prandtl number /3-convection forced by thermal 
boundary fluxes or internal heating is an area for future study. 

With respect to Jovian dynamics, a problem with the model is that it 
generates only one zonal shear zone, and multiple jets in y have not yet 
been found (though see fig. 6b for some wrinkles). Soward (1977) investi- 
gated the dynamics of small but finite-amplitude convection in an inter- 
nally heated fluid sphere and concluded that meridional variations should 
appear on scales of order 0 -2 / 9 as /3 -*■ «, which only becomes substantially 
less than one at the extreme values of 0 in our calculations. It is not 
clear how relevant this weakly nonlinear analysis is for the turbulent 
regime we are studying. For this situation Rhines (1975) has proposed that 
2-D turbulence on a /3-plane will cascade energy upscale into zonally symme- 
tric modes with y-wavenumbers k*^ « r/3*/2U*rms]'/ 2 , or equivalently, into 
non-dimensional cell sizes 1^ ■ */Dk ^ = 2 3 / 4 ir K' '!*/&'/* , Using the data 
from table 1 it can be seen that 1^ is of order 1 except for those runs with 
0 - in where the convection is only weakly turbulent. In order to study 
the possibility of multiple jet solutions we need to consider simulations 
at even higher R a , perhaps as large as 10 10 . Although the possibility of 
attaining flow statistics over a wider range of parameters is appealing, 
such simulations will require much higher resolution and very tiny time 
steps . 

The two-dimensional /3-convection model has allowed us to study flow 
evolution and equilibration at higher Rayleigh numbers and Taylor numbers 
(« /3 2 ), and over much longer time scales, than is possible with a global 
eddy-resolving 3-D convection code. However, as computer parallelism and 
speed increase, it would be interesting to look for similar turbulent vac- 
illation states in three dimensional compressible models of stellar and 
planetary atmospheres. 
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FIGURE CAPTIONS ' 

Fig. 1. Geometry of the equatorial channel model. The channel has width D 
and height L. Gravity is assumed perpendicular to the rotation axis. 

Fig. 2. Properties of linear instabilities, (a) Critical Rayleigh number 
vs. zonal wavenumber (lowest cross-stream mode) for values of 0 
shown, (b) Growth rate Re(s) for various /3 at R a = 10® . 

Fig. 3. Typical wavenumber spectra of nonlinear solutions. 0 = 2800. (a) 

Symmetric wavenumber 1, R a = 42,500. (b) Asymmetric wavenumber 1, R a 
= 75000, with comparable energies in low and moderate wavenumbers, 
(c) Asymmetric type 2A(6) with dominant wavenumber 6, R a = 57500. (d) 

Asymmetric type 8A, R a = 40000. 

Fig. 4. State diagram for 0 = 2800, tracking various wavenumber and symme- 
try states vs. Rayleigh number. Each horizontal line represents a 
solution branch with F = steady amplitude travelling waves (fixed 
point in energy), L = limit cycle, L2 = period-doubled cycle, T = 
two-frequency torus, T2 = period-doubled 2-torus, C = chaos. 
Single-tipped arrows denote changes upon increasing Rayleigh number 
using the last timestep of the previous solution as the initial condi- 
tion. Double-tipped arrows denote evolution upon small amplitude 
perturbations (=1 percent) with no symmetry. 

Fig. 5. Results along the Type 8A branch. 0 - 2800. (a) R a = 57,500. (b) R a 
= 75,500. (c) R a = 76,000. (d) R a = 125,000. Vertical columns from 
left to right show area-averaged zonal kinetic energy vs. its rate of 
change, Poincare' sections of the phase plots, time series of zonal 
kinetic energy, and frequency spectra of the times series. 

Fig. 6. Eddy (m ?£ 0) stream function and zonal velocity vs. time for a) case 
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corresponding to Fig. 5a. b) corresponding to Fig. 5d. For this and 
all shaded pictures the gray scale spans the range of negative or min- 
imum (black), to positive or maximum (white) values. 

Fig. 7. Results along the Type IS (unstable) branch. 0 = 2800. (a) R a = 

43.000. (b) R a = 43,250. (c) R a = 43,600. (d) R a = 43,750. The struc- 

ture of the diagram is the same as Fig. 5. 

Fig. 8. Total stream function and horizontal wavenumber spectra for the 
multiple states at R a = 60,000, 0 = 2800. (a) The chaotic type 1A(6) 

mode, (b) the steady amplitude travelling wave 1A(7) mode. 

Fig. 9. Times series of eddy (m *0) streamf unction at R a = 10* . (a) 0 = 0. 
(b) 0 = 6000 (ordinary convection branch), (c) 0 = 6000 (/3-convection 
branch), (d) 0 = 8000. (e) 0 = 20,000. (f) 0 = 28,000. Coloration spans 
the range of negative (black) to zero (red) to positive (white) values. 

Fig. 10. Dominant wavenumbers at R a = 10*. The crosses indicate spectral 
peaks, while the squares show the wavenumber of maximum linear 
growth rate. 

Fig. 11. Variation of space averaged zonal kinetic energy and Nusselt 
number with time for the large zonal jet cases from figure 9 at R a = 
10*. The Nusselt number is 8T/dy at y = 0. 

Fig. 12. Pulsating zonal jets and turbulent convection at R a = 10 7 , 0 = 

75.000. Total stream function and non-conduct ive temperature are 
illustrated with the color scale as in figure 9. 

Fig. 13. Time series of total mean kinetic energy and the pointwise zonal 
velocity at y * 0.6667 for the case shown in figure 12. The Nusselt 

number here is ^ at y = 0. The portions between the vertical dashed 

lines are illustrated in figure 12. 

Fig. 14. Ratio of time averaged mean flow kinetic energy to time averaged 
eddy kinetic energy over the whole channel for the R a values shown. 
The horizontal portions, past the /3-cutoff (the short dashed curve), 
show the zero-lines for these offset plots. The long-dash curve in- 
dicates the transition from ordinary convection with very weak mean 
flows to those with large zonal jets. The narrow window of hyster- 
esis in the transition is illustrated for R a = 10*. 

Fig. 15. Qualitative regime diagram showing locations of pure conduction 
('Jstable"), /?-convection with large zonal jets, and "ordinary" sta- 
tionary or travelling Rayleigh-Benard convection with small or zero 
(if 0 = 0) zonal flows. The asterisk points locate parameters where 
the largest zonal jets are found, subject to our fairly coarse parame- 
ter space sampling. The time-average cross-channel velocity differ- 
ences are indicated. 
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